home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / timidsrc.zip / filter.c < prev    next >
C/C++ Source or Header  |  1996-05-20  |  6KB  |  204 lines
  1. /*
  2.  
  3.     TiMidity -- Experimental MIDI to WAVE converter
  4.     Copyright (C) 1995 Tuukka Toivonen <toivonen@clinet.fi>
  5.  
  6.     This program is free software; you can redistribute it and/or modify
  7.     it under the terms of the GNU General Public License as published by
  8.     the Free Software Foundation; either version 2 of the License, or
  9.     (at your option) any later version.
  10.  
  11.     This program is distributed in the hope that it will be useful,
  12.     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13.     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14.     GNU General Public License for more details.
  15.  
  16.     You should have received a copy of the GNU General Public License
  17.     along with this program; if not, write to the Free Software
  18.     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  19.  
  20.    filter.c: written by Vincent Pagel ( pagel@loria.fr ) 
  21.  
  22.    implements fir antialiasing filter : should help when setting sample
  23.    rates as low as 8Khz.
  24.  
  25.    April 95
  26.       - first draft
  27.  
  28.    22/5/95
  29.       - modify "filter" so that it simulate leading and trailing 0 in the buffer
  30.    */
  31.  
  32. #include <stdio.h>
  33. #include <string.h>
  34. #include <math.h>
  35. #include <stdlib.h>
  36. #include "config.h"
  37. #include "common.h"
  38. #include "controls.h"
  39. #include "instrum.h"
  40. #include "filter.h"
  41.  
  42. /*  bessel  function   */
  43. static float ino(float x)
  44. {
  45.     float y, de, e, sde;
  46.     int i;
  47.     
  48.     y = x / 2;
  49.     e = 1.0;
  50.     de = 1.0;
  51.     i = 1;
  52.     do {
  53.     de = de * y / (float) i;
  54.     sde = de * de;
  55.     e += sde;
  56.     } while (!( (e * 1.0e-08 - sde > 0) || (i++ > 25) ));
  57.     return(e);
  58. }    
  59.  
  60. /* Kaiser Window (symetric) */
  61. static void kaiser(float *w,int n,float beta)
  62. {
  63.     float xind, xi;
  64.     int i;
  65.     
  66.     xind = (2*n - 1) * (2*n - 1);
  67.     for (i =0; i<n ; i++) 
  68.     {
  69.         xi = i + 0.5;
  70.         w[i] = ino((float)(beta * sqrt((double)(1. - 4 * xi * xi / xind))))
  71.         / ino((float)beta);
  72.     }
  73. }
  74.  
  75. /*
  76.  * fir coef in g, cuttoff frequency in fc
  77.  */
  78. static void designfir(float *g , float fc)
  79. {
  80.     int i;
  81.     float xi, omega, att, beta ;
  82.     float w[ORDER2];
  83.     
  84.     for (i =0; i < ORDER2 ;i++) 
  85.     {
  86.         xi = (float) i + 0.5;
  87.         omega = PI * xi;
  88.         g[i] = sin( (double) omega * fc) / omega;
  89.     }
  90.     
  91.     att = 40.; /* attenuation  in  db */
  92.     beta = (float) exp(log((double)0.58417 * (att - 20.96)) * 0.4) + 0.07886 
  93.     * (att - 20.96);
  94.     kaiser( w, ORDER2, beta);
  95.     
  96.     /* Matrix product */
  97.     for (i =0; i < ORDER2 ; i++)
  98.     g[i] = g[i] * w[i];
  99. }
  100.  
  101. /*
  102.  * FIR filtering -> apply the filter given by coef[] to the data buffer
  103.  * Note that we simulate leading and trailing 0 at the border of the 
  104.  * data buffer
  105.  */
  106. static void filter(sample_t *result,sample_t *data, int32 length,float coef[])
  107. {
  108.     int32 sample,i,sample_window;
  109.     int16 peak = 0;
  110.     float sum;
  111.  
  112.     /* Simulate leading 0 at the begining of the buffer */
  113.      for (sample = 0; sample < ORDER2 ; sample++ )
  114.     {
  115.         sum = 0.0;
  116.         sample_window= sample - ORDER2;
  117.        
  118.         for (i = 0; i < ORDER ;i++) 
  119.         sum += coef[i] *
  120.             ((sample_window<0)? 0.0 : data[sample_window++]) ;
  121.         
  122.         /* Saturation ??? */
  123.         if (sum> 32767.) { sum=32767.; peak++; }
  124.         if (sum< -32768.) { sum=-32768; peak++; }
  125.         result[sample] = (sample_t) sum;
  126.     }
  127.  
  128.     /* The core of the buffer  */
  129.     for (sample = ORDER2; sample < length - ORDER + ORDER2 ; sample++ )
  130.     {
  131.         sum = 0.0;
  132.         sample_window= sample - ORDER2;
  133.  
  134.         for (i = 0; i < ORDER ;i++) 
  135.         sum += data[sample_window++] * coef[i];
  136.         
  137.         /* Saturation ??? */
  138.         if (sum> 32767.) { sum=32767.; peak++; }
  139.         if (sum< -32768.) { sum=-32768; peak++; }
  140.         result[sample] = (sample_t) sum;
  141.     }
  142.     
  143.     /* Simulate 0 at the end of the buffer */
  144.     for (sample = length - ORDER + ORDER2; sample < length ; sample++ )
  145.     {
  146.         sum = 0.0;
  147.         sample_window= sample - ORDER2;
  148.         
  149.         for (i = 0; i < ORDER ;i++) 
  150.         sum += coef[i] *
  151.             ((sample_window>=length)? 0.0 : data[sample_window++]) ;
  152.         
  153.         /* Saturation ??? */
  154.         if (sum> 32767.) { sum=32767.; peak++; }
  155.         if (sum< -32768.) { sum=-32768; peak++; }
  156.         result[sample] = (sample_t) sum;
  157.     }
  158.  
  159.     if (peak)
  160.     ctl->cmsg(CMSG_ERROR, VERB_NORMAL, 
  161.           "Saturation %2.3f %%.", 100.0*peak/ (float) length);
  162. }
  163.  
  164. /***********************************************************************/
  165. /* Prevent aliasing by filtering any freq above the output_rate        */
  166. /*                                                                     */
  167. /* I don't worry about looping point -> they will remain soft if they  */
  168. /* were already                                                        */
  169. /***********************************************************************/
  170. void antialiasing(Sample *sp, int32 output_rate )
  171. {
  172.     sample_t *temp;
  173.     int i;
  174.     float fir_symetric[ORDER];
  175.     float fir_coef[ORDER2];
  176.     float freq_cut;  /* cutoff frequency [0..1.0] FREQ_CUT/SAMP_FREQ*/
  177.  
  178.  
  179.     ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: Fsample=%iKHz",
  180.           sp->sample_rate);
  181.  
  182.     /* No oversampling  */
  183.     if (output_rate>=sp->sample_rate)
  184.     return;
  185.     
  186.     freq_cut= (float) output_rate / (float) sp->sample_rate;
  187.     ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: cutoff=%f%%",
  188.           freq_cut*100.);
  189.  
  190.     designfir(fir_coef,freq_cut);
  191.     
  192.     /* Make the filter symetric */
  193.     for (i = 0 ; i<ORDER2 ;i++) 
  194.     fir_symetric[ORDER-1 - i] = fir_symetric[i] = fir_coef[ORDER2-1 - i];
  195.     
  196.     /* We apply the filter we have designed on a copy of the patch */
  197.     temp = safe_malloc(sp->data_length);
  198.     memcpy(temp,sp->data,sp->data_length);
  199.     
  200.     filter(sp->data,temp,sp->data_length/sizeof(sample_t),fir_symetric);
  201.     
  202.     free(temp);
  203. }
  204.